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This  report  summarizes  the  analyses  of  deflagration  to  detonation 
transition  [DDT)  occurring  in  a  packed  bed  of  granular,  high  energy  solid 
propellant.  A  reactive  two-phase  flow  model  of  this  phenomena  is  solved 
by  utilizing  a  lax-Wendroff  finite  differencing  technique.  A  brief 
overview  of  the  well  known  shock  jump  conditions  for  one-dimensional, 
one-phase  flow  with  heat  addition  is  reported,  and  a  similar  analysis 
for  one-dimensional ,  two -phase  reactive  flow  is  discussed.  Improvements 
made  in  the  gas  phase  nonideal  equation  of  state,  gas  permeability,  and 
numerical  integration  t  chniques  allow  for  the  prediction  of  a  transition 
to  a  steady  detonation  from  initiation  by  deflagration. 

Analyses  are  presented  that  clearly  indicate  the  effect  of  the 
propellant  physical  and  chemical  parameters  on  the  predicted  run-up 
length  to  detonation.  Predictions  of  this  run-up  length  to  detonation 
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CHAPTER  ONE 


DEF  L AGRAT I ON -^0 - DETONAT I ON  TRANS I T ION 

1 . 1  Introduction 

This  report  summarizes  the  analysis  associated  with  the  accelerating 
deflagration  wave  in  a  porous  medium  of  reactive  solid  propellant.  The 
phenomenon  of  DDT  (.def lagration-to-detonation  transition)  in  solid  pro¬ 
pellants,  e  pecially  solid  propellants  burning  in  rocket  motor  environ¬ 
ments,  is  not  usually  considered  a  hazard.  However,  it.  may  be  that  un¬ 
der  certain  situations,  for  example  a  grain  structure  failure,  the  solid 
motor  may  crack  and  form  regions  of  granular  or  porous  propellant.  When 
flame  from  the  surface  deflagrating  propellant  reaches  this  seam  of  porous 
material  it  will  accelerate  into  this  medium  and  be  supported  by  convective 
heat  transfer  from  the  burnt  gas  into  the  unignited  porous  region.  If  in 
addition  to  this  the  product  gases  are  confined  to  a  finite  volume,  the 
accelerating  deflagration  could  transit  into  a  detonation.  Propellants 
exhibiting  a  high  chemical  energy  per  unit  mass  and  capable  of  rapidly 
generating  gases  through  their  burning  rate  are  more  likely  to  experience 
this  type  of  deflagration  to  detonation  transition. 

Analysis  of  the  flow.*  in  such  an  unsteady  two-phase  mixture  is  a  com¬ 
plicated  exercise.  Work  has  been  underway  at  the  University  of  Illinois 
since  197S  to  develop  a  reactive  hydrodynamic  code  in  which  the  combus¬ 
tion  of  porous  propellants  can  be  modeled  in  such  a  way  as  to  predict  the 
behavior  of  a  convectively  driven  flame  in  a  confined  situation.  Details 
of  these  modeling  exercises  are  found  in  References  1-4.  In  an  evolutionary 
mariner  this  work  included  the  formulation  of  the  two-phase  flow  conserva¬ 
tion  equations,  first  assuming  that  the  mixture  was  a  continuum,  and  at  a 
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later  date  treating  each  phase  a:,  a  separate  continuum  irrespective  oz 
the  mixture  properties.  The  most  recent  analysis  of  the  unsteady  two-phase 
flows  associated  with  DDT  is  documented  in  the  paper  by  Hoffman  and  Krier 
[3].  This  reference,  therefore,  represents  the  smarting  point  for  the 
work  that  is  presented  here. 

The  reader,  after  reviewing  the  above  noted  reference,  will  understand 
that  the  modeling  of  this  transient  phenomena  utilizes  a  number  of  impor¬ 
tant  constitutive  relations,  which  form  closure  of  the  conservation  equa¬ 
tions.  For  example,  one  must  have  information  on  the  burning  rate  of  the 
material  that  is  a  function  of  the  surrounding  pressure  and  temperature. 

One  must  also  have  laws  for  the  dynamic  gas  permeability  and  the  subse¬ 
quent  heat  transfer  rates  of  the  hot  gases  as  they  are  forced  into  the 
unignited  porous  material.  In  addition,  relations  which  represent  the 
resistance  to  compaction  of  the  solid  matrix  must  be  included.  Equations 
of  state,  not  only  for  the  high  pressure  in  the  product  gases,  but  also 
for  the  solid  itself,  must  be  supplied  (as  shown  in  the  work  by  Hoffman 
and  Krier).  The  assumption  of  an  incompressible  solid,  although  provid¬ 
ing  some  reasonable  answers  s  far  as  the  deflagration  speed,  is  not  an 
accurate  indicator  of  the  peak  pressures  that  are  possible  during  the 
accelerating  deflagration  mode.  Since  these  pressures  are  precursor  to 
the  final  detonation  solutions  that  would  be  expected,  it  is  clear  that 
a  compressible  solid  must  be  modeled. 

In  summary,  the  analysis  of  the  DDT  problem  requires  the  solution  of 
the  conservation  equations  in  both  phases  and  the  necessary  constitutive 
relations.  The  conservation  equations  form  a  system  of  nonlinear  hyper¬ 
bolic  equations,  which  require  numerical  finite  differencing  schemes. 

This  report  will  summarize  several  of  these  integration  schemes  and  will 
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evaluate  which  are  more  useful  for  this  reactive  flow  problem.  There  are, 
of  course,  many  numerical  techniques  available  to  handle  hyperbolic  par¬ 
tial  differential  equations  which  allow  for  solutions  of  shock  waves. 
Obviously,  not  all  of  these  have  been  treated  in  this  work. 

1 . 2  Previous  Results  on  DDT  Modeling 

A  review  of  References  2-4  indicates  that  a  steady  state  detonation 
solution  was  not  a  predicted  result.  As  will  be  shown  in  Sections  1.3 
and  2.4,  for  propellant  chemical  energy,  the  initial  solids  loading  con¬ 
sidered  and  the  material  parameters  (e.g.  y) ,  a  steady  state  detonation 
(CJ)  would  propagate  at  speeds  of  the  order  of  5-8  nun/ysec,  with  a  deto¬ 
nation  pressure  of  the  order  of  12-20  GPa .  Although  a  fairly  rapid  flame 
front,  often  approaching  a  steady  state  speed  of  2  mm/ysec  was  typical  of 
the  solutions  presented  in  the  work  by  Gokhale  and  Krier  [2] ,  Kezerle  and 
Krier  [4],  and  Hoffman  and  Krier  [3],  the  associated  peak  pressures  were 
never  of  significant  manner  to  suggest  that  these  final  flame  spreading 
rates  were  characteristic  of  an  actual  detonation,  since  fronts  ranging 
from  j  to  2  mm/ysec  cannot  support  a  detonation  (see  following  section) . 
However,  the  peak  pressures  previously  reported,  which  ranged  between 
1/2  to  5  GPa,  were  consistent  with  the  velocities. 

Continued  work,  which  is  reported  here,  indicates  that  it  is  now 
possible  to  obtain  a  detonation  solution,  but  in  order  to  do  so  a  num¬ 
ber  of  important  modifications  and  corrections  are  required.  To  under¬ 
stand  what  these  changes  are  and  why  they  are  necessary,  it  is  appro¬ 
priate  at  this  point  to  first  provide  a  review  on  the  topic  of  detona¬ 
tion.  Although  the  discussion  that  follows  is  documented  in  various 
manners  in  several  textbooks,  including  such  a  review  here  will  allow 
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for  a  clearer  understanding  of  the  steady  state  solution  that  is  termed 
"detonation",  as  well  as  define  the  flow  properties  of  the  DDT  analysis. 

1 . 3  Jump  Conditions  for  Reactive  Flow 

In  order  to  understand  the  jump  conditions  across  a  detonation  wave 
in  two-phase  (sol?d-gas)  reactive  flow,  one  should  first  look  at  the  solu¬ 
tion  of  the  one-phase  (all  gas]  flow  problem.  Pickett  and  Davis  [5],  anu 
Strehlow  [6]  -'evelop  these  jump  equations  in  more  detail  in  their  texts. 
This  section  will  give  the  reader  a  brief  review  of  their  work  and  high¬ 
light  some  of  the  key  assumptions  made  in  developing  the  Hugoniot  and 
Rayleigh  line  equations. 

In  Figure  1  a  detonation  wave  with  velocity  D  is  moving  through  a 
constant  area  duct  into  the  unignited  or  "cold"  end  which  is  at  rest. 

For  this  analysis  all  chem  cal  reactions  are  assumed  to  take  place  in  a 
narrow  reaction  zone.  The  products  of  combustion  are  shown  moving  with 
velocity  u  in  the  same  direction  as  D.  In  this  figure,  the  subscript  A 
on  pressure,  density,  and  velocity  designates  the  unignited  or  reactants 
side,  and  B  indicates  the  product  side  of  the  combustion  zone. 

To  better  understand  the  jump  conditions  across  the  shock  in  Figure 
1,  a  new  frame  of  reference  may  be  helpful.  From  a  coordinate  system 
attached  to  the  moving  detonation  wave,  an  observer  located  at  the  origin 
would  see  the  reactants  moving  with  velocity  u^  =  D  into  the  stationary 
reaction  zone  and  the  products  of  combustion  exiting  with  velocity 
Ug  =>  (D-u)  (see  Figure  2). 

Conservation  of  ass  and  momentum  through  the  flow  area  shown  in 
Figure  2  gives  respectively 

PAU  "  PBUB  =  PB(D“U) 


(1) 
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Moving  Reaction  Zone7 


Fig.  1.  Detonation  Wave  Moving  into  Stationary  Reactants. 
(Reference  Frame  Fixed) . 
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In  this  model  all  viscous  effects  are  confined  within  the  bounds  of  the 
shock  discontinuity.  The  expression  for  conservation  of  momentum  (Eq.  2) 
can  be  simplified  by  making  use  of  Equation  1  to  obtain 
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PA  '  V” 


(5) 


Later,  it.  will  be  assumed  that  PA  is  nagligible  with  respect  to  P^,  but 
for  the  time  PA  will  be  retained  in  the  momentum  equation. 

Elimination  of  u,  the  reaction  products  velocity  expressed  in  the 
fixed  coordinate  system,  from  Eqs .  1  and  3  gives  the  well  known  Rayleigh 
or  Michelson  line. 

p/d"  -  (PD  -  Pa)/(va  -  vB)  *  0  (4) 


Figure  3  illustrates  a  Rayleigh  line  on  a  P-v  diagram  for  one  specific 

2  2 

value  of  D.  The  line  has  slope  -pA  D  and  passes  through  the  point  C?A»  VA) • 

For  the  case  where  D  -*■  the  Rayleigh  line  approaches  vertical  and  for 

D  -*■  0,  the  Rayleigh  line  is  horizontal.  In  order  to  satisfy  the  conserva¬ 
tion  of  mass  and  momentum,  the  final  state  fP  .  v  1  must  also  lie  nn  the 

. '  '  '  “  B'  '  B'  •” 

Rayleigh  line  as  shown  in  Figure  3. 

Simultaneously  solving  Equations  1  and  3  by  eliminating  D  instead  of 
u,  the  result  would  be  an  expression  for  u,  the  gaseous  product  velocity 
in  the  fixed  reference  frame. 


uZ  *  (PB  -  PA)  (VA  -  V  (*) 

Figure  4  shows  several  curves  of  constant  u  plotted  on  a  P-v  diagram. 
Solving  Equations  1  and  2  for  =  (D-u)  one  would  obtain 
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Fig,  3.  Illustration  of  Rayleigh  lir.e. 
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(5a) 


Equation  5a  will  be  important  later  in  showing  that  for  a  CJ  detonation 
the  product  velocity  is  sonic  with  respect  to  the  detonation  front. 

The  third  conservation  equation,  the  conservation  of  energy  across 
the  reaction  zone,  can  be  written 


E.  „  +  P.v.  +  4  D2  *  E.  „  +  P.vD  +  i  uD2  =  E.  „ 

int.  A  A  2  intQ  B  B  2  B  mtn 

A  B  B 


*  PBVB  *  I  (D'U)2 


(6) 


In  Equation  6,  represents  the  specific  internal  energy  at  each  re¬ 
spective  location  and  v  is  the  specific  volume  v  =  1/p.  Since  Cv>  the 
specific  heat  at  a  constant  volume,  is  assumed  constant  throughout  the 
process,  the  internal  energy  at  location  B  is  simply 


E.  =  C  T_ 
v  8 


(7a) 


Likewise,  the  internal  energy  at  location  A  is  similar  with  the  inclusion 
of  the  chemical  energy  yet  to  be  released. 


E.  =  C  T  +  E 
int.  v  A  CHF.M 
A 


(7b) 


By  convention  a  positive  value  of  ^  is  endothermic  and  a  negative  value 


ovorhormi < 


Elimination  of  u  and  D  from  the  energy  equation  (Eq.  6)  is  obtained  by 
substituting  in  the  mass  and  momentum  equations  (I'qs.  1  and  3).  The  re¬ 
sult  is  a  relation  between  t.  P,  and  v  in  both  reactant  and  prcduct 

int’  r 

states  known  as  the  Hugoniot  equation.  It  is  expressed  as 


tint  2  (PB 
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V  = 


(8) 


If  one  assumes  that  is  a  constant  and  the  reaction  products  obey  an  ideal 
equation  of  state 
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Pv  =  RT 

then  the  internal  energy  expressions  (Eqs.  7a  and  7b)  take  the  form 


(9) 
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(10b) 


where  y  is  the  ideal  gas  limit  of  Substitution  of  Equations  10a 

and  10b  into  Equation  8  gives  a  new  form  of  the  Hugoniot  equation. 
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On  a  P-v  diagram  Equation  11  plots  as  a  hyperbola  passing  through  (P  ,  v^) 
for  the  case  where  E^^  =  0  and  displaced  from  this  point  for  increasing 
values  of  E^^  (see  Fig.  5).  From  the  solution  of  the  conservation  equa¬ 
tions,  the  Hugoniot  relation  states  that  flow  coming  in  at  pressure  P  =  P 

A 

and  specific  volume  v  =  VA  with  chemical  energy  E^^  =  ECHEM  ,  must  have 

a  final  state  (P  ,  v  )  which  lies  on  the  Hugoniot.  As  discussed 

o  D  LrihM^ 

earlier,  the  Rayleigh  line  connecting  the  initial  state  (P  ,  v  )  with  the 

A  JK 

final  state  (Pg,  v^)  dictates  a  unique  solution  for  the  detonation  velocity, 
D,  for  that  process.  Therefore,  the  steady  s  ate  flow  process  between 
states  A  and  B  shown  in  Figure  2  can  be  uniquely  modeled  by  the  Rayleigh 
line  and  Hugoniot  curve  on  Figure  6,  knowing  the  initial  state  (P, ,  v  ), 
the  final  stat<  (Pg,  Vg)  and  the  chemical  energy  of  the  re  .ctant. 

There  are  certain  areas  on  the  P-v  diagram  where  a  final  solution  to  the 
flow  process  is  .impossible.  For  instance,  the  pressure  and  specific 


Hugomot 
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volume  cannot  both  increase  over  the  shock  discontinuity . 

Figure  6  shows  a  Kugoniot  curve  and  two  Rayleigh  lines,  and  P7, 
on  the  same  P-v  diagram.  For  Rayleigh  line  labeled  R2  there  are  two  solu¬ 
tions  labeled  S  for  strong  and  W  for  weak.  At  the  strong  solution  the 
flow  downstream  of  the  shock  (uD  in  Fig.  2)  is  subsonic  with  respect  to 

D 

the  stationary  shock,  and  a  pressure  disturbance  initiated,  traveling  at 
the  local  speed  of  sound  in  the  fluid,  will  propagate  back  and  overtake 
the  shock  front.  Looking  at  the  same  strong  solution  but  in  the  fixed  re¬ 
ference  frame  (moving  detonation  wave) ,  a  pressure  fluctuation  behind  the 
detonation  front  (i.e.,  moving  the  rear  wall  in  a  closed  system)  will  over¬ 
take  the  front  and  the  front  will  adjust  itself  to  that  change.  For  the 
weak  solution  the  flow  behind  the  front  (ug  in  Fig.  2)  moves  away  from 
die  stationary  front  at  a  velocity  greater  than  the  local  sound  velocity. 
Therefore,  a  pressure  disturbance  in  this  case  will  not  be  felt  by  the 
front.  Both  the  strong  and  weak  detonations  will  be  discussed  in  more 
detail  in  the  next  section. 

In  the  case  of  the  Rayleigh  line  labeled  R^  and  the  Hugoniot  curve, 
there  is  only  one  unique  solution  to  the  conservation  equations.  At  this 
point  the  two  lines  are  tangent  and  the  final  st.  te  is  identified  as  the 
CJ  (Chapman- Jouguet)  point.  For  this  case,  the  flow  behind  the  front  at 
the  CJ  condition  is  moving  at  the  local  speed  of  sound  with  respect  to 
the  front  (i.e.,  ufi  =  a  )  . 

In  order  to  prove  that  the  flow  is  sonic  with  respect  to  the  deto¬ 
nation  front, one  should  first  look  at  the  thermodynamic  definition  of 
sound  speed . 


(12) 


In  Equation  12  the  subscript  s  indicates  the  derivative  is  evaluated 
at  constant  entropy.  Therefore,  for  a  detonation  the  local  sound  velocity 
a_  the  CJ  point  is  simply  obtained  by  evaluating  the  slope  of  an  isen- 
trepe  through  that  point  on  a  P-v  plot. 

Since  the  slope  of  the  Hugoniot  and  Rayleigh  lines  are  equal  at  the 
CJ  point,  it  is  useful  to  solve  both  derivatives  and  equate  them.  For 
the  Rayleigh  line 

<f>  ■  -  °Ad2  ■  -  (PB  -  V-'CvA  -  V  <13> 

K 

and  for  the  Hugoniot  curve 

<f>H  •  2‘f>n  '  tvA  -  V  ‘  <PB  -  PJ  '  ‘"A  -  V  t») 

Equating  the  right  hand  sides  of  Eqs.  13  and  14  you  obtain 

"  -  P  U5) 

HCJ 

From  the  thermodynamic  relation 


dE  =■  Tds  -  Pdv 


(16) 


you  obtain 


(17) 


by  evaluating  Equation  16  on  an  isentrope.  Equations  15  and  17  imply 
that  at  the  CJ  point  an  isentrope  is  tangent  to  both  Rayleigh  and  Hugoniot 
lines.  From  this 
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Therefore,  at  the  CJ  point  the  local  sound  speed  can  be  evaluated  by  using 
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liquations  12  and  18  to  obtain 
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This  is  exactly  the  same  as  Equation  5a,  the  square  of  the  product  ve¬ 
locity  for  a  stationary  detonation.  Therefore,  the  product  velocity  for 
a  CJ  detonation  is  sonic  with  respect  to  the  detonation  front. 

From  Equation  18,  one  can  define  gamma  to  be  <the  negative  logarith¬ 
mic  slope  of  the  isentrope  through  the  CJ  point. 


y  =- (3ln  P/31n  v)  = 

The  analysis  of  detonation  processes  presented  in  this  section  as¬ 
sumes  y  is  not  only  a  constant  value  but  is  the  same  in  both  product  and 
reactant  sides  of  the  detonation  wave.  For  most  gases  Y  is  about  1.2  [7], 
In  the  case  of  a  sc  lid  (explosive)  detonating,  Fickett  and  Davis  [5]  use 
a  value  of  y  =  3.0  and  assumes  the  ideal  state  equation  holds  for  the 
product  gases.  From  this  assumption  and  estimate  of  y,  one  could  directly 
apply  the  analysis  formulated  in  this  section  to  the  solid  detonation. 

However,  for  two-phase  reactive  flow  the  mixture  momentum  and  mix¬ 
ture  energy  differ  from  Eqs.  3  and  4  by  interactive  stress  and  stress 
work  terms.  Mso,  the  gaseous  products  cannoc  be  assumed  ideal  when  ex¬ 
perimental  work  shows  detonation  pressures  of  the  order  15-20  GPa  (2.17-2.9 
Mpsi)  [8] .  Because  of  the  high  detonation  pressure  a  nonideal  equation 
of  state  v, ; s  implemented.  This  will  be  discussed  in  more  detail  in  Chapter 
2.  Also,  the  jump  conditions  for  a  mixture  with  interactive  terms  is 
outlined  in  Appendix  A.  It  is  similar  to  the  maverial  presented  in  this 
section  with  the  exception  of  the  stress  and  stress  work  terms. 
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1.4  Rear  Boundary  Condition  for  Steady  State  Detonations 


Reference  5  shows  that  the  velocity,  D,  at  which  a  detonation 
propagates  into  the  reactants  is  dependent  on  the  velocity  at  which  the 
rear  wall,  or  boundary,  travels  with  respect  to  the  CJ  product  velocity 
u  .  The  product  velocity  was  discussed  in  the  previous  section  on  the 
Rayleigh  line  and  Figure  4  shows  several  curves  of  constant  product  ve¬ 
locity.  The  CJ  product  velocity  would  be  the  curve  intersecting  the  CJ 
point  in  Figure  6.  Three  different  possibilities  for  the  movement  of 
the  rear  wall  with  respect  to  the  CJ  product  velocity  exist  and  will  be 
discussed  in  this  section. 

The  first  case  to  consider  is  where  the  rear  wall  velocity  is  greater 
than  the  product  velocity  at  the  CJ  condition.  On  the  P-v  diagram  (see 
Fig.  6)  this  would  correspond  to  the  solution  labeled  S  for  "strong". 
Here,  the  flow  field  following  the  detonation  is  constant  with  the  pro¬ 
duct  velocity  u  =  uw>  Examining  Figure  6,  one  can  see  that  the  final 
state  pressure,  P_,  is  also  greater  than  that  of  the  corresponding  CJ  con- 
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dition.  Because  the  flow  behind  the  detonation  front  is  subsonic  with 
respect  to  the  front,  any  perturbation  in  the  rear-  wall  will  propagate 
forward  to  the  pressure  front.  Figure  7a  shows  a  pressure  profile  for  an 
overdriven  detonation  with  the  CJ  pressure  labeled  on  the  figure. 

The  eff'  cts  on  the  pressure  by  reducing  the  rear  wall  velocity  can 
be  seen  in  Figure  7b.  Here,  a  rarefaction  from  the  decreased  wall  pres¬ 
sure  is  shown  propagating  forward  to  the  detonation  front.  An  actual 
overdriven  or  strong  detonation  must  have  an  external  force  driving  the 
rear  wall  at  a  velocity  greater  than  the  CJ  product  velocity.  Fickett 
and  Davis  [5]  give  as  an  example  for  an  overdriven  detonation  the  case 
where  another  detonation,  stronger  than  the  one  being  studied,  is  driving 
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the  rear  wall. 

A  second  possibility  would  be  the  case  where  the  wall  velocity  is 
exactly  equal  to  the  product  velocity  at  the  CJ  condition.  Here,  the 
detonation  pressure  is  the  CJ  value  and  the  flow  following  the  front  is 
sonic  with  respect  to  the  front.  Figure  8  illustrates  case  2. 

The  third,  and  most  important  case  for  our  study  of  DDT  initiated 
by  flame  in  a  closed  tube,  is  when  the  wail  velocity  is  less  than  the  CJ  ve- 
1  city.  For  this  case,  the  front  still  propagates  at  the  minimum  detona¬ 
tion  velocity  allowable  by  the  conservation  equations,  the  CJ  detonation 
velocity,!^,  and  the  combustion  products  leave  at  velocity  ,  sonic  with 
respect  to  the  front.  Since  the  wall  velocity  is  less  than  the  product  ve¬ 
locity,  a  rarefaction  from  the  rear  of  the  combustion  zone  to  the  wall  is 
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This  is  an  iscH tropic  expansion  characterized  by  the  ratio  of 


specific  heat,  y,  assumed  for  the  products  (see  discussion  in  previous  sec¬ 
tion).  Figure  9  illustrates  the  case  where  the  wall  velocity  is  less  than 
the  product  velocity.  As  an  example  of  the  isentropicity  cf  the  flow  be¬ 
hind  the  front,  the  pressure  ratio  between  the  final  reaction  state  (CJ) 
and  the  location  labeled  k  in  Figure  9  is  obtained  from  the  expansion 
wave 


l +  O^r  )  M 


T2y/(Y-1) 

kJ 


where  M.  is  the  Mach  number  of  the  flow  at  location  k. 
k 
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1 . 5  Steady  and  Unsteady  Deflagration 

Along  with  the  three  detonation  solutions  of  the  Rayleigh  line  and 
Hugoniot  curves  (i.e,  CJ,  strong  solution,  and  weak  solution)  discussed 
in  the  previous  section,  there  are  ether  steady  state  solutions  possible. 
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Figure  10  illustrates  a  solution  to  the  conservation  equations  on  the  lower, 
or  deflagration,  branch  of  the  Hugoniot  curve.  For  this  solution  the 
flame  propagates  through  the  unburnt  reactants  at  a  low  Mach  number  (M«l) 
and  as  can  be  seen  on  the  figure,  the  deflagration  is  characterized  by  a 
drop  in  pressure  and  rise  in  specific  volume  across  the  reaction  zone.  Re¬ 
call,  for  the  detonation  solution  there  was  a  rise  in  pressure  and  d„op  in 
specific  volume.  An  example  of  a  steady  state  deflagration  would  be  the 
case  of  a  flame  burning  on  a  bunscn  burner.  Here,  the  flame  is  stationary 
and  the  reactant  gas  is  flowing  through  the  flame  sheet  at  a  low  velocity. 

When  discussing  DDT  in  a  packed  porous  bed,  one  must  clarify  exactly 
the  fluid  dynamics  of  the  deflagration  branch.  At  time  t=Q,  one  end  of 
the  packed  bed  is  ignited  and  generating  gas  from  the  propellant  surface. 

As  time  progresses,  it  is  the  hot  gases  confined  in  the  region  behind  the 
ignition  front  which  in  a  sense  drives  the  ignition  front,  through  the  bed. 

That  is,  unlike  che  pressure  drop  across  the  steady  state  deflagra¬ 
tion  solution,  there  is  a  pressure  rise  across  the  ignition  front.  The 
increase  in  pressure  is  a  result  of  the  gas  generation  in  the  region  be¬ 
hind  the  ignition  front  being  confined  to  a  finite  volume.  Also,  in  the 
analysis  made  of  the  steady  state  deflagration,  all  reaction  was  assumed 
to  take  place  across  an  infinitely  thin  re^rcion  zone.  I*  is  obvious 
that  for  the  two-phase  deflagration  problem  this  is  not  the  case.  Be¬ 
cause  of  the  finite  zone  of  reaction,  the  problem  of  the  deflagration 
flame  is  at  most  quasisteady  and  it  may  be  tha*  the  nonsteady  terms  in 
the  conservation  equations  should  not  be  neglected.  Because  of  this 
inherent  unsteadiness  of  the  deflagration  phase  of  the  DDT  problem,  the 
reader  should  not  confuse  it  with  the  steady  state  deflagration  solut'on 
shown  in  figure  10.  Later  it  will  be  shown  that  if  transition  from 


22 


deflagration  to  detonation  does  occur,  the  reaction  zone  must  collapse 
and  a  detonation  solution  can  be  analyzed  as  being  steady  and  by  the 
jump  conditions  being  properly  satisfied. 

1 . 6  DDT  in  Two-Phase  Reactive  Flow  (Experimental  Data) 

As  discussed  in  the  previous  section,  the  transition  from  deflagra¬ 
tion  to  detonation  in  a  porous  reactive  medium  is  an  unsteady  process. 

Hot  gases  generated  from  the  propellant  surfaces  are  driven  forward  into 
the  unburnt  solid  matrix  by  the  pressure  gradient  developed  at  the  igni¬ 
tion  front.  This  phenomena  is  not  found  when  a  nonporous  solid  deto¬ 
nates,  since  only  pressure  disturbances  can  be  propagated  ahead  of  the 
ignition  front.  It  is  this  convective  heat  transfer  to  the  unignited 
propellant  and  the  extended  deflagration  reaction  zone  which  makes  DDT 
in  a  porous  bed  unique  when  comparing  it  to  DDT  in  either  an  all  gas  or 
an  all  solid  regime. 

Much  experimental  work  on  DDT  in  porous  material  has  beer,  carried 
out  in  the  past  decade.  Bernecker  and  Price  [8-10],  have  published  the 
most  recent  results  on  DDT  in  a  series  of  three  papers.  Other  experi¬ 
mental  studies  prior  to  these  include  the  work  of  Griffiths  and  Grocock 
[11]  and  Taylor  [12]. 

The  work  presented  in  Reference  9  by  Bernecker  and  Price  is  a  study 
on  DDT  in  RDX  ( cyclotrimethylenetrinitramine) ,  a  shock  sensitive  high 
energy  explose.  In  their  experiments  the  RDX  was  packed  into  a  thick 
walled  tube  having  an  inside  diameter  of  16mm  and  being  approximately 
300mm  (12  in.)  in  length  (see  Fig.  lia  [8]).  Both  ends  of  the  column  were 
closed  and  ionization  probes  were  located  throughout  tne  bed  to  trace  the 
ignition  front  locus.  The  RDX  was  packed  in  an  inert  wax  mixture  and 
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had  a  mean  particle  size  of  200ym  in  diameter.  Figure  lib  illustrates  the 
DDT  mechanism  on  a  distance-time  plot  for  the  91/9  RDX/wax  granular  charge 
[9].  4.s  shown  in  the  figure,  the  convective  flame  front  obser/ed  in  their 

work  had  no  acceleration  up  to  the  detonation  transition. 

Their  experimental  work  showed  this  convective  ignition  front  traveled 
at  subsonic  velocities  [0. 3-0. 9mm/ u sec)  for  most  cases  when  the  initial 
porosities,  <t>o>  ranged  between  values  of  0.1  to  0.3.  The  lower  data  points 
on  the  plot  represent  a  post  convective  compressional  wave  in  the  burning 
region  which  overtakes  the  ignition  front.  The  length  of  porous  propel¬ 
lant  it  takes  for  the  compressional  wave  to  overtake  the  ignition  front 
and  transit  into  a  detonation  is  i^j,  the  run-up  length. 


1 . 7  Topics  to  be  Addressed 

The  work  nresented  in  this  chanter  reviewed  those  features  which  dis- 

A  A 

tinguish  DDT  in  a  porous  reactive  medium  from  DDT  in  other  media.  Since 
the  work  that  follows  in  s^usequent  chapters  will  attempt  to  model  this 
phenomena,  it  was  important  that  the  properties  of  an  actual  detonation 
be  understood(Sections  3-5) . 

It  was  also  pointed  out  that  the  DDT  process  being  studied  cannot  re¬ 
present  the  transition  between  the  lower  CJ  point  and  the  uppor  CJ  point. 
That  is,  because  the  flame  propagates  from  a  closed  end  and  the  deflagra¬ 
tion  reaction  zone  is  never  of  infinitesim  1  thickness,  one  cannot  apply 
steady  state  jump  conditions  for  the  deflagration  phase. 

The  chapter  that  follows  reviews  the  past  DDT  modeling  efforts  and 
inidcates  how  this  work  builds  for  the  present  effort.  Numerical  solu¬ 
tions  to  the  flow  equations  are  carried  out  am  solutions  are  presented 
in  Chapter  4.  These  results  will  show  that  the  detonation  (steady  state 
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supersonic  wave)  is  a  solution  for  a  certain  class  of  problems.  Such  re¬ 
sults  have  yet  to  be  presented  by  others. 


I 
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CHAPTER  TWO 

THE  MODEL  AND  ANALYSIS 


2. 1  Introduction 

The  analysis  that  follows  attempts  to  model  the  situation  in  which 
a  bed  of  tightly  packed  granular  propellant  is  ignited  at  one  end.  Both 
ends  of  the  packed  bed  are  considered  closed,  thus,  modeling  the  problem 
discussed  in  Section  1.4  where  the  velocity  of  the  rear  wall  was  less 
than  the  CJ  product  velocity.  Recall  for  this  case  that  after  DDT,  the 
detonation  propagates  through  the  unignited  region  at  the  CJ  detonation 
velocity  and  detonation  pressure  and  is  followed  by  an  isentropic  expan¬ 
sion  of  the  gases  to  the  pressure  at  the  stationary  wail.  The  gas  sur¬ 
rounding  the  particles  at  the  initial  time  is  considered  inert,  and  to 
be  at  atmospheric  pressure.  It  is  also  assumed  that  the  inert  gas  will 
fully  mix  with  the  gases  being  generated  from  combustion  of  the  propel¬ 
lant  as  time  progresses.  Figure  12  is  a  schematic  of  the  propellant  bed. 

For  numerical  simplicity  the  propellant  particles  are  assumed  to  be 
unisized  spheres.  Particles  of  interest  range  in  diameter  50pm  <  d^  <  200pm. 
To  treat  mult  -sized  particles,  one  would  require  _N  independent  equations 
of  mass,  momentum  and  energy  for  the  solid  phase,  whore  N  is  the  number 
of  initially  diff crent-uized  particles.  A  solids  loading  of  74%  is  the 
tightest  possible  for  uni  sized  spheres,  obtainable  by  arranging  the  spheres 
in  a  face  centered  cubic.  However,  assuming  granular  deformation  occurs 
under  high  stress  loads,  as  discussed  by  Kuo  [13 j  ,  greater  solids  loadings 
may  be  predicted  without  error.  Obviously,  the  spherical  geometry  must 
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he  altered  for  this  to  occur. 

As  the  small  fraction  of  propellant  particles  ignited  at  time  t=0 
burn,  hot  gases  are  generated  as  a  function  of  the  pressure-dependent 
burning  rate  law  and  surface-to-volume  ratio  (3/r  )  of  the  spheres.  These 
hot  gases  generated  are  convected  forward  through  the  lattice  of  unburned 
propellant  and  flow  gradients  develop,  as  dictated  by  the  solution  of  the 
conservation  equations  and  the  necessary  constitutive  relations. 

Heat  transfer  from  the  hot  gases  to  the  unignited  propellant,  par¬ 
ticles,  dependent  on  the  velocity  of  the  gas  relative  to  the  particles 
and  several  gas  properties  (i.e.  viscosity,  thermal  conductivity),  trans¬ 
ports  energy  from  the  gas  to  the  solid  phase.  Subsequent  ignition  of  par¬ 
ticles  further  down  the  bed  is  assumed  to  occur  when  a  critical  solid  phase 
internal  energy  is  reached  [2],  This  energy  can  be  expressed  as  a  critical 
increase  in  solid  phase  temperature,  T  ,  since  the  specific  heat  of  the 
solid  is  assumed  to  be  known. 

As  time  progresses  the  gas  pressure  behind  the  ignition  front  increases 
due  to  the  confinement  of  the  gases  from  the  closed  rear  boundary  and  the 
pressure-dependent  rate  of  mass  generation  in  the  gas  phase.  Under  certain 
conditions*  the  pressure  gradient  can  develop  into  a  shock  front  which 
overtakes  the  ignition  frcnt  propagating  through  the  bed.  When  this  occurs 
the  ignition  front  experiences  the  transition  from  deflagration  to  detonation 
discussed  in  Section  1.6. 

At  the  transition  point  the  ignited  region  (zone  of  gas  generation) 
narrows  in  width  and  is  followed  by  a  region  of  all  gas  where  the  propellant 


These  conditions  will  depend  upon  the  solid  chemical  energy,  granula¬ 
tion  (size  and  loading),  burning  rate,  ignition  energy,  etc. 
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particles  are  completely  burned  out  (see  Fig.  13).  The  thickness  of  the 
reaction  zone  is  a  function  of  the  initial  particle  size  and  solids  load¬ 
ing,  and  thicknesses  approaching  1mm  may  be  possible  as  r^  +  0  and  4>0  -*■  0. 


2 .  2  Assumptions 

In  order  to  numerically  model  DDT  in  two-phase  flow  while  retaining 
the  physics  of  the  problem,  several  key  assumptions  had  to  be  made.  These 
assumptions  are  similar  to  those  made  by  previous  investigators  (Ref.  2 


and  4  ) . 


(1)  Both  the  solid  and  gas  phases  are  independently  treated  as  con¬ 
tinuums  requiring  their  own  conservation  relations. 

(2)  Each  phase  interacts  with  the  other.  This  is  modeled  by  the 
mass,  momentum  and  energy  interaction  terms  in  the  conserva¬ 
tion  equations. 

f 3)  All  propellant  particles  are  unisized  spheres. 

[4)  Ignition  of  a  propellant  particle  is  obtained  when  a  critical 
energy,  expressed  as  a  particle  temperature,  is  transferred  to 
the  solid. 


1 I ant  particles  are  initially  surrounded  by  an  inert 


gas  at  temperature,  T 

*o 

(6)  During  combustion  of  the  propellant,  the  gaseous  products  mix 
with  the  inert  gas  described  in  assumption  (5) . 

(7)  Both  ends  of  the  bid  are  closed  allowing  no  gases  to  escape. 

(8)  The  specific  heats  at  constant  volume,  ,  for  both  phases  are 
constant . 


(.9)  When  the  solid  phase,  at  a  given  x- location  in  the  bed,  displays 
a  porosity  <p  >0.95  and  a  particle  volume  one-tenth  of  its  initial 


volume,  it  is  burned  out  and  no  longer  generating  gas.  Results 
show  this  phenomena  to  proceed  smoothly  from  the  left  boundary 
and  thus  not  leaving  any  'holes'  in  the  continuum.  This  assump¬ 
tion  was  necessary  to  prevent  a  singularity  from  arising  as  r^  -*■  0 
and  T  -*•  00 . 

(10)  All  the  product  gases  obey  an  assumed  nonideal  equation 
of  state. 

(11)  The  solid  particles  are  compressible,  without  heating  up,  obey¬ 
ing  a  modified  Tait  equation  of  state. 

(12)  Cnee  ignited,  the  particles  are  assumed  to  burn  on  the  outer 
surface  only,  at  a  known  pressure-dependent  rate  law. 

(13)  At  some  initial  time,  a  "narrow”  region  at  one  end  is  ignited, 
burning  at  the  low  pressure  prescribed. 

2 . 3  Governing  Equations 

Numerical  modeling  of  the  two-phase  reactive  flow  process  in  DDT  in¬ 
volves  the  conservation  of  mass,  momentum  and  energy  per  unit  volume  in 
both  gas  and  solid  phases.  This  is  a  system  of  six  conservation  condi  - 

+■  i  rtnc  uVi  i  r>Vi  4-  /-\  r*»n  n  cat  n  ■£"  nnnl  i  naer  Vm  fv*  av-Uri  1  i  r>  -I  o  1  n  1 

***v**j  *  vim  w  -Jvw  w*.  uv/aaw  put  umoi  uiibiai  CVJUO~ 

tions  coupled  by  the  interphase  mass,  momentum  and  heat  transfer  terms.  The 
conservation  equations  for  two-phase  reactive  flow  have  been  developed 
previously  and  References  3  and  4  will  provide  the  reader  with  the  defi¬ 
nitions,  assumptions  and  expressions  for  the  six  following  field-balanced 
conservation  eauations.  The  e  are: 

Gas  Continuity 
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Particle  Continuity 
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Particle  Energy 


a  CM)  iiP2y  *  O-WVp1 


at 


dx 


f£p _ H] 

1  CHEM  2  J 


+  V  u  •*  Q 
P 


(27) 


Here,  the  relations  for  the  total  internal  energy  in  each  pnase  are 

E<t  5  Cvg  Tg  *  I  “*  a"d  V  VP‘‘*  "p 

The  subscripts  g  and  p  denote  gas  and  particle  ■  espectively.  In  Equations 
22  to  27,  the  phase  densities,  and  p^»  are  defined  as 


P!  =  pg* 


and 


P2  =  (l-4>)pp 


The  porosity,  <J>,  is  defined  as  the  ratio  of  the  instantaneous  gas  volume 

Hence,  the  solids  fraction  is  (!-$)• 


to  the  mixture  volume. 
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In  addition  to  the  six  conservation  equations,  three  constitutive 

relations  are  needed  in  order  to  solve  for  the  nine  unknown  variables; 

p,p,u,u,T,T,<t>,P  and  P  .  These  relations  include  state  equa- 
g  P  g  P  g  P  g  P 

tions  for  both  gas  and  solid  states  and  a  stress-resistance  relation  for 
Pp.  Appendix  B  gives  a  complete  listing  of  the  relations  used. 

2.4  Improvements 

Since  the  work  reported  in  Reference  3,  certain  "improvements"  in 
the  modeling  effort  have  allowed  solutions  which  may  be  considered  to  be 
actual  detonations.  These  improvements  are  discussed  in  some  detail  later 
in  the  text,  but  basically  include: 

1.  Implementation  of  the  necessary  gas  phase  (nonideal)  equation 
of  state,  to  insure  that  at  the  CJ  (Chapman-Jouget)  cc.  ditions, 
the  isentrope  provides  for  a  "gamma  law"  suitable  at  the  hydro- 
dynamic  CJ  state.  (Appendix C  presents  a  review  of  such  an 
equation  of  state.) 

2.  Implementation  of  a  new  gas-particle  friction  coefficient,  a: 
developed  by  Wilcox  and  Krier  [14],  for  flows  at  the  high  Reynolds 
numbers  encountered  in  the  developing  DDT  flows.  Previously  such 
coefficients  were  based  on  data  only  available  fos  moderate  Reynolds 
number  ranges. 

3.  Utilization  of  a  modified  numerical  integration  scheme,  which 
allowed  for  a  reduction  in  the  grid  spacing  (and  hence  reduc¬ 
tion  of  the  time  increment)  without  the  usual  penalty  of  exces¬ 
sive  computation  costs  ana  numerical  instability  that  often  follows 
when  the  total  number  of  integrations  is  significantly  increased. 
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2 . 5  Equations  of  State:  Constraint  Due  to  the  Detonation  State 

As  mentioned  in  Chapter  1,  a  nonideal  equation  of  state  must  be  uti¬ 
lized  for  the  product  gases.  The  analysis  presented  here  uses  a  nonideal 
equation  of  state  for  hard  spheres  developed  by  S.  J.  Jacobs  [15].  Pre¬ 
vious  to  this  study,  a  covolume-type  state  equation  with  data  made  avail¬ 
able  by  Cook  [16]  was  used  [see  discussion  in  Ref.  2) . 

The  hard  sphere  equation  of  state  takes  the  form 


1.0  +  bp  +  c  (bp)  +  .  .  . 


(28a) 


where  the  constants  b  and  c  are  determined  by  the  value  of  the  gamma  law 
coefficient,  y,  for  the  product  gases  as  discussed  in  Appendix  B.  As 
stated,  y  is  the  negative  logarithmic  slope  of  the  lsentro^e  tangent 
to  both  t’ne  Rayleigh  and  Ilugonrot  lines  "t.  the  CJ  point  m  the  detona¬ 
tion  state.  That  is,  the  slope  of  the  isentrope  that  the  product  gases 
expand  along  in  the  product  state.  The  reader  may  refer  to  Appendix  C 
for  the  complete  solution  of  the  constants  b  and  c. 

Values  for  y  range  from  two  to  three  for  detonating  high  density 
explo  ives  [15].  For  the  baseline  case  considered  in  this  study,  a 
value  of  y  =  2.05  was  selected  and  the  corresponding  nonideal  equation  of 


state  is: 


~r~  =  1.0  *  2.5  p  -  0.50  p  2  (28b) 

RT  g  g 

When  the  above  coefficients  (2.5  and  0.50)  were  altered  to  treat  a  case 
for  y  =  3.0,  excessively  high  gas  and  particle  temperatures  were  predicted 
as  one  mignt  expect.  In  addition,  during  the  numerical  integration  severe 
oscillations  in  ;as  and  particle  temperatures  occurred  in  most  cases  when 
y  >  3.  Since  one  is  always  constrained  by  the  numerical  integration  schemes 


36 


that  are  employed  to  handle  very  severe  gradients  in  the  flow,  it  is  under¬ 
standable  that  previous  efforts  in  DDT  modeling  [2,  3],  which  without 
knowing  were  utilizing  high  y  values  for  the  product  gases,  almost  always 
ran  into  numerical  failure. 

Evaluating  the  covolume  state  equation  previously  used  in  Ref.  3,  an 
approximate  value  of  y  =  3.6  is  calculated.  This  extreme  value  for  y  may 
be  one  important  reason  why  the  calculations,  as  reportel  in  Ref.  2-4, 
were  unable  to  handle  the  high  pressures  associated  with  steady  state 
detonations . 

In  the  solid  phase  the  particles  obey  a  modified  Tait  equation  allow¬ 
ing  for  compression  of  the  granules.  This  is  written  as: 


0 


K 


1/3 


+  1 


(29) 


*  r  o  *-  o  -» 

wnere  Kq  is  the  bulk  modulus.  Reference  3  discusses  the  Tait  equation  in 
more  detail.  A  typical  value  is  Kq  =  1.38  GPa  (2.0  x  IQ5  psi) . 


2 . 6  GJ  (Detonation)  State 

In  the  text  by  Fickett  ard  Davis  [5]  equations  are  presented  for  es¬ 
timating  detonation  velocity,  D  .  ,  and  the  detonation  nressure.  p_.  .  for 
a  steady  state  detonation.  These  are  given  for  a  single  phase  explosive: 


and 


[CJ  =  2Cr‘1-)Pp0  echem 


dcj=  echem 


(30) 


(31) 


Here,  p  is  the  initial  solid  material  density  and  E  is  the  chemical 
p0  CHEM 

energy  liberated  by  burning  the  solid.  Since  the  problem  being  considered 
is  two-phase  (solid-gas),  Eqs.  30  and  31  must  be  modified  to  account  for 
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this  by  converting  p  ,  the  initial  solid  density,  to  p~  ,  the  initial 

P0  o 

solid  phase  density  3  P„  (1-40  •  Foi  an  initial  solid  density  p^  =  1994 


o  *  o 

=0.30,  a  chemical  energy  E 
2.05,  Eqs.  30  and  31  give  respectively: 


kg/m  ,  an  initial  porosity  4>rt  =0.30,  a  chemical  energy  E 
and  assuming  y 


=5.48  MJ/kg, 


and 


CJ 


CJ 


16.06  GPa 


5.92  -HU 
qsec 


These  equations  were  developed  from  the  jump  equations  for  one-phase  flow 
where  the  equation  of  state  for  the  product  gases  was  assumed  ideal.  Be¬ 
cause  of  this,  these  equations  should  only  be  used  to  get  a  good  esti¬ 
mate  of  the  detonation  pressure  and  velocity. 

In  addition  to  Equations  30  and  31  another  method  was  presented  by 
Fickett  and  Davis  [5]  for  estimates  of  detonation  pressure  and  detona¬ 
tion  velocity.  This  is  hamlet' s  Short  Method  and  was  developed  by  hamlet 
and  Jacobs  [17],  in  the  CJ  state  they  are: 


0.5 


PCJ  = 

P 


1000 


(32) 


and 


0.25 


°CJ-  Atl*Boo  )n"W  •31-62  l3i) 

P 

where  the  constants  A  =  2.23  (m-kg-s  ^(mole  MJ)  B  =  0.0013 

(m3/kg)  and  K  =  0.762  (Nm^  kg  *(mole  MJ)  ,  Again,  p  is  the  initial 

3  P 

solid  density  in  kg/m  and  N  the  reciprocal  of  the  gaseous  molecular 

weight  in  mole/kg.  Equations  32  and  33  give  detonation  pressure  and  ve¬ 
locity  for  a  detonating  solid  explosive  going  to  all  gas  in  the  product 

state.  To  modify  Equations  32  and  33  for  the  two-phase  problem, p  ,  the 

D 

solid  density,  is  again  converted  to  the  solid  phase  density  P0  =  (l-4>  )p 

Pc 
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It  should  be  noted  that  Equation  33  expresses  detonation  velocity  as 
a  function  of  propellant  density,  while  Equation  31  shows  it  to  be  inde¬ 
pendent  of  density.  Like  Equations  30  and  31,  Equations  32  and  33  should 
only  be  taken  as  approximations. 

Using  the  same  input  as  above.  Equations  32  and  33  give  respectively 


CJ 


CJ 


»  16.52  GPa 


=  7.83 


mm 

usee 


2 . 7  Gas  Permeability 

One  of  the  key  constitutive  relations  required  in  the  analysis  is 
the  gas-particle  (interphase)  viscous  force,  which  governs  hot  gas  pene¬ 
tration  into  the  unignited  region  of  the  granular  material.  As  presented 
in  Appendix  A: 

U 

V  =  — ^  (u  -  u  )  f  (34) 

4r  2  g  P  Pg 

P 

where  f  is  the  drag  coefficient.  Until  recently,  the  packed  bed  corre¬ 
lations  by  Ergun  or  Kuo  and  Nydegger  (as  reviewed  in  Ref.  14)  were  uti¬ 
lized  for  f  .  Thus,  the  modeling  efforts  presented  in  Ref.  3  and  4  used 
the  expression  of  Kuo  and  Nydegger  [18] : 

i2  .81 


f  .  {276  +  5  (-—-)  } 


Pg 


l-d 


(35) 


Equation  35  was  developed  for  460  <  Re  <  14,600.  Here,  Re  is  the  appro¬ 
priate  Reynolds  number,  defined  as: 


Re  -  [(<t>ug)pg  2rp]/Pg  (36) 

Based  on  experiments  at  both  high  gas  velocities  and  high  Reynolds  numbers, 
Wilcox  and  Krier  [14]  developed  the  correlation: 
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f  =  5.06  x  10S  r  Re/u2 
Pg  P  8 


(37) 


where  u  is  given  in  m/sect  r  in  meters,  and  the  constants  5.06  x  10 
g  s  P 

must  have  units  of  m/s2.  This  expression  is  not  valid  for  either  very 

low  gas  velocities,  since  Equation  37  would  give  f  ■*  00  as  -*•  0,  or 

for  very  high  gas  velocities*,  since  f  -*■  0  as  u  -*•  ®.  The  equation 

3  5 

was  found  to  be  fairly  accurate  for  10  <  Re  <  (2  x  10  )and  15  m/s<  u  <  150  m/s. 

g 

While  a  straight  forward  comparison  is  not  easy,  the  difference  between  the 

value  for  f  as  predicted  by  Equation  35  versus  that  by  Equation  37  can 

4 

be  seen  in  the  example.  At  Re  =  10  and  $  =  0.4,  Equation  35  gives 

f  =  53,600  and  Equation  37  gives  a  value  of  f  =  5448.  To  utilize 
Pg  H  Pg 

Equation  37  one  must  specify  the  average  gas  velocity,  u  ,  and  the  par- 

g 

tide  radius,  ro,  that  were  used  to  obtain  the  Reynolds  number.  A  kine- 

“6  2 

matic  viscosity  y  / p  =  1.8  x  10  m  /sec,  a  particle  radius  r  =  1.0mm, 
g  g  P 

and  an  average  gas  velocity  u^  =  30.5  m/sec  were  used  to  obtain  a  Reynolds 
~  4 

number  Re  =  10  .  From  this  particular  example,  the  Wilcox/ Krier  correla¬ 
tion  (Eq.  37)  allows  about  10  times  the  permeability,  i.e.,  1/10  the  viscous 
drag  force  as  correlated  by  the  Kuo/Nydegger  relation. 

9  'll  Hi'S  nnc  tn  flow  nr-nr  <»<; «;  1  uadi  no  fn  DOT  di  rn<;  «tA<'  in  Tnllnw- 

- -  —  - - -  r - - - - ©  - *  - - - -  -  - -  - 

ing  chapter,  clearly  indicate  that  sufficient  gas  permeability  is  necessary 

to  allow  for  a  detonation  transition. 


•Data  from  Ref.  14  was  limited  to  u  < 


300  m/ sec . 
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CHAPTER  TmEE 

NUMERICAL  INTEGRATION 


3. 1  Finite  Difference  Mesh 

To  solve  the  Eulerian  formulated  system  of  conservation  equations 
discussed  in  Chapter  ?  with  the  constitutive  relations  listed  in  Appendix 
B,  the  length  of  the  bed  being  integrated  over  is  divided  into  I  segments, 
each  a  constant  Ax  in  width  (i.e.,  x.  =  jAx;  j  =  1,  2,  3.,..  I).  The  value 
required  for  Ax  will  be  discussed  later  in  the  chapter. 

At  t  =  0  values  of  the  nine  independent  variables;  p  ,  p  ,  u  ,  u  , 

*  8  P  8  P 

T  ,  T  ,  P  ,  P  and  $  are  initialized  at  each  ith  x-location  in  the  erid. 
g  P  g  P 

Before  incrementing  the  primary  variables  (i.e.  mass,  momentum  and  energy) 
to  the  future  time,  t  -  tQ  +  At,  the  auxiliary  variables  in  the  equations 
(e.g.,drag,  gas  generation,  heat  transfer)  must  be  computed  at  the  present 
time,  t  =  t  .  The  nine  equations  are  then  solved  at  the  incremental  time, 
t  =  t  +  At, by  a  modified  Lax-Wendroff  finite  differencing  scheme.  This 
method  along  v  .th  another  tested  are  presented  in  Appendix  D.  The  second 
method,  developed  by  Rubin  and  Berstein  [19],  was  implemented  for  several 
test  cases  and  proved  to  give  results  similar  to  that  of  the  Lax-Wendroff 
scheme.  However,  to  keep  all  results  consistent,  the  Lax-Wendroff  scheme 
was  cho’-en  ?.s  the  integrator. 

The  time  Increment,  At,  over  which  the  equations  are  solved  is  cal¬ 
culated  by  the  Cour&nt,  Fredrich;; ,  Levy  stability  criteria  [20],  for  hy¬ 
perbolic  equations, 

XAx 

s  (Fi^n 


At 


(38) 


1 

1 
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In  Equation  38,  the  term  c  is  the  mixture  sound  speed  and  |u|  is  the  maxi¬ 
mum  gas  velocity  in  the  bed.  Also.  X  is  a  stability  constant,  less  than 
unity.  For  most  cases  X  =  0.5  was  used.  This  smaller  time  increment 
was  necessary  in  order  to  integrate  the  equations  ’.'hen  large  gradients 
developed  in  the  flow. 

3 . 2  Initial  end  Boundary  Conditions 

To  initialize  the  problem,  the  bed  is  assumed  to  be  quiesent,  i.e., 
at  a  constant  gas  temperature  and  constant  gas  pressure  throughout  the 
length  of  the  bed.  The  spherical  propellant  particles  are  typically  fix¬ 
ed  at  a  constant  solids  loadings,  although  a  variable  initial  porosity 
can  be  treated.  Then  to  initiate  the  flow,  the  propellant  at  the  first 
few  grid  points  is  assumed  ignited  and  generating  gas.  To  be  consistent 
with  the  fiat  initial  pressure  profile,  all  gas  and  particle  velocities 
at  time  t  =  0  are  set  equal  to  zero. 

In  the  past  an  exponential  pressure  profile  was  constructed  at  the 
initial  time  in  order  to  speed  up  the  deflagration  process.  This  was 
inaccurate,  however,  in  that  it  was  unsupported  (i.e.  ,  no  fluid  motion 
was  associated  with  the  pressure  gradient).  Nevertheless,  the  computa¬ 
tion  time  necessary  for  the  pr  'ssure  gradient  to  develop  on  its  own  was 
too  great  and  the  initial  pressure  gradient  was  implemented  to  speed  up 
the  proces-'  Since  then  the  computer  code  was  modified  to  include  an  in¬ 
tegration  process  which  cuts  the  computation  .ime  drastically  and  allows 
for  the  use  of  a  more  realistic  flat  initial  pressure  profile  (P  ~  P  ) . 
Table  1  is  a  summary  of  typical  input  data  for  the  cases  studied. 

3 .  3  Modification  to  the  1  n  t  egration  Scheme 

To  solve  on  a  digital  nputer  the  flow  which  is  to  represent  a  DDT 
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Table  1 

TYPICAL  INPUT  DATA 


_ _ PARAMETER _ 

Burning  Rate  Index 

Burning  Rate  Proportionality  Constant 

Initial  Bed  Porosity 
Particle  Diameter 
Bed  Length 
Grid  Spacing 
Chemical  Energy 

Gas  Specific  Heat 


VALUE 


0.8  <  n  <  1.2 

b  =  0.001  —  [~1 

sec  lPSIJ 


0.25  <  4>  <  0.50 

—  o  — 

50um  dp  <_  SOOum 
L  =  25. 0  cm 


Ax  =  1 . 27  mm 

4  —  <  If  <  7  Md 

kg  CHEM  7  kg 

If  T 

L0kA^cv^L9 


KJ 

kg^k 


Initial  Bed  Temperature 
Ignition  Temperature  (Bulk) 

Ignition  Energy 


T  294  °K 

g 

T.  =  303  °K 
ign 

E .  => 

ign 


9.0 


KJ 

H 
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phenomena,  requires  the  repetitive  integrations  of  many  equations.  All 
of  the  conservation  equations  and  constitutive  relations  are  solved  indi¬ 
vidually  at  each  grid  point  i\i  the  entire  bed,  for  each  time  increment. 

A  study  of  several  test  cases  showed,  for  instance,  that  when  the  ig¬ 
nition  front  was  at  a  given  x- location,  there  was  little  activity  several 
grid  points  ahead  of  it.  That  is,  the  gas  and  particles  were  at  velocities 
close  to  zero  and  all  transport  coefficients  were  negligable  in  the  region 
not  far  in  front  of  the  ignition  front.  This  phenomena  occurred  for  both 
the  detonation  state  where  the  front  was  propagating  between  5-8mm/usec, 
and  the  deflagration  state  where  the  propagation  velocity  was  much  less. 
Because  of  this,  an  addition  was  made  to  the  computer  code  which  allows 
the  code  to  only  integrate  the  active  region  and  bypass  the  inactive  zone. 

Ahead  of  the  ignition  zone,  the  computer  code  locates  the  nearest 
point  to  the  zone  where  there  is  no  significant  particle  or  gas  movement. 
This  point  is  then  designated  as  the  new  front  boundary  of  the  integration 
region  for  that  particular  integration  step.  When  the  pressure  front 
builds  the  ignition  front  moves  rather  rapidly  through  the  bed  (5-10-™^.-) 
and  a  new  integration  boundary  must  be  located  after  each  time  increment. 

Figures  14  and  15  show  s  comparison  of  a  case  run  without  the  moving 
integration  zone  (Fig.  14)  and  one  with  (Fig.  15).  To  assure  a  correct 
solution  the  actual  front  integration  boundary  is  extended  a  few  grid 
points  beyend  the  location  calculated  by  the  code.  The  addition  of  this 
new  logic  to  the  current  code  reduced  computation  ti  e  by  at  least  one 
half,  for  the  same  time  increment  and  grid  space. 

3 . 4  Artificial  Smoothing 

Inherent  to  the  solution  of  the  system  of  interdependent  conservation 


equations  is  a  numerical  instability.  A  small  perturbation  can  in  some 
cases  amplify  with  each  time  increment  and  eventually  destroy  the  numeri¬ 
cal  solution.  This  phenomena  can  start  at  the  first  time  increment  and 
in  ten  to  twenty  integrations  the  oscillations  can  be  so  large  that  the 
solution  becomes  unstable  and  terminates.  In  order  to  smooth  out  these 
oscillations  before  they  amplify,  an  artificial  smoothing  routine  must 
be  incorporated  into  the  code. 

From  experience  in  integrating  the  two-phase  flow  equations  (Eqs. 
22-27)  with  significant,  nonlinear  source-sink  terms,  the  problem  of  nu¬ 
merical  instability  occurs  often  enough  to  warrant  artifical  smoothing. 

An  extensive  study  by  the  author  and  other  investigators  [2-4] ,  shows  the 
final  solution  to  be  independent  of  any  artificial  smoothing  used- 

The  analyzation  of  numerous  test  cases  has  shown,  suprisingly.  that 
the  stable  solution  did  not  require  smoothing  of  all  the  variables.  It 
is  obvious  that  the  following  results  have  an  inherent  dependency  on 

smoothing  .  Extreme  care  has  been  taken  to  minimize  these  effects  on  the 
qualitative  trends  and  quantitative  results  predicted.  Nevertheless,  it 
should  be  obvious  that  smoothing  techniques  can  supply  variability  in  the 
predicted  parameters  which  are  reflections  of  the  scheme  utilized,  and  not 


necessarily  of  the  conservation  equations. 

The  particular  artificial  smoothing  technique  found  to  be  useful  is 
a  three  point  averaging  technique  where  the  variable  v  determines  what 


4? 


whore 


For  example,  if  v  =  0.1  then  W  retains  80%  of  its  original  value  and  is 
influenced  by  10%  of  the  U  values  on  each  side.  This  is  shown  graphically 
in  Figure  16. 

For  all  cases  studied,  a  minimum  value  of  V  was  used  that  would  pro¬ 
vide  for  a  stable  solution.  This  required  repeating  each  test  case  nume¬ 
rous  times,  lowering  v  ea?h  time  until  the  solution  went  unstable. 

3. 5  Grid  Spacing 

It  is  obvious  that  in  order  to  minimize  computation  time,  a  maximum 
value  of  Ax  which  gives  a  stable  solution  should  be  used.  Based  on  the 
calculations  carried  out  in  References  2-4  a  grid  space  of  Ax  =■  1.27  mm 
was  shown  to  be  the  largest  Ax  that  would  provide  a  solution  to  a  flow 
problem  which  exhibited  rather  steep  gradients. 
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CHAPTER  FOUR 


RESULTS  COMPUTED 


4 . 1  Introduction 

This  chapter  will  present  the  calculations  made  on  the  possibility 
of  DDT  occurring  in  packed  beds  of  high  energy,  granulated,  unisized  pro¬ 
pellants  or  explosives.  It  is  obvious  that  there  are  an  infinite  number 
of  loading  combinations  possible  for  the  DDT  study.  Fortunately,  the  work 
of  Hoffman  and  Krier  [3]  ana  Krier  and  Gokhale  [2],  in  which  conditions 
of  rapid  convective  flame  spreading  have  been  calculated,  is  available 
and  can  be  used  as  a  starting  point.  It  should  be  noted  that  none  of 
the  calculations  made  in  References  2  or  3  predicted  DDT.  As  pointed  out 
in  the  previous  chapters,  there  are  probable  reasons  why  this  has  not  been 
accomplished  and  it  is  expected  that  certain  improvements  and  modifica¬ 
tions  will  now  allow  for  the  calculation  of  the  steady  state  detonation 
solution. 

The  study  made  of  a  DDT  potential  attempts  to  model  conditions  simi¬ 
lar  to  the  test  conditions  of  the  Bernecker  and  Price  work  [9]  fi.e.,  a 
long  column  of  granulated  material  in  a  closed  pipe  ignited  by  an  ener¬ 
getic  ignition  material  at  one  end] .  In  order  to  model  the  flow  tran¬ 
sient,  one  must  assure  that  the  lengtn  of  the  bed  exceeds  £  ,  the  run-up 

length  to  detonation.  Since  the  experimental  work  of  Bernecker  and  Price 
[9]  indicated  that  a  10  in  (25  cm]  bed  was  sufficient  for  most  of  their 
experiments  where  DDT  occurred,  this  length  was  selected  as  the  longest 
bed  length  to  be  considered. 
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Since  unisized  spheres  are  being  treated,  the  initial  porosity  can 
be  no  less  than  $  =  0.26,  although  randomly  packed  unisized  spheres 

generally  give  a  high  porosity,  about  =  0.40.  Therefore,  for  this 
study  0.26  <_  <  0.40.  The  initial  particle  radii  studied  were  also  in 

the  same  range  as  those  considered  by  References  2  and  3.  Results  in 
these  studies  showed  that  particles  must  be  less  than  one  millimeter  in 
diameter  in  order  to  generate  sufficient  gases  for  the  rapid  flame  spread¬ 
ing  phenomena. 

The  chemical  energy  of  the  material  considered  is  in  the  range  of 
explosives  or  high  energy  propellants  of  the  nitramine  family  (i.e.,  HMX, 

RDX) .  Thus,  the  chemical  energies  studied  were  always  larger  than  1000 
cal/g  (4.15  MJ/kg) .  Other  parameters  one  must  consider  in  the  DDT  studies 
are  the  burning  rate  properties  of  the  propellant.  Again,  the  values  used 
in  References  2  and  3,  which  attempted  to  model  the  burning  rate  of  an 
HMX  solid  propellant,  were  utilized.  However,  the  burning  rate  index, 
n,  is  a  parameter  which  is  explicitly  studied. 

In  this  analysis,  the  deflagration  will  be  initiated  by  assuming  at 
time,  t  =  0,  that  except  for  a  small  portion  of  the  bed  at  one  closed  end, 
the  bed  is  quiescent  and  unreacting.  As  has  been  mentioned  in  Chapter  3, 
a  closed  end  situation  is  considered  and  hence  at  the  two  end  points  (  x  =  0, 
x  =  L)  it  is  assumed  that  all  flow  gradients  are  zero  and  that  the  veloci¬ 
ties  of  the  particles  and  gas  must  also  be  zero. 

The  following  section  begins  with  a  solution  that  indicates,  for  the 
first  time,  a  steady  state  detonation  can  be  predicted.  This  result  is 
then  compared  to  conditions  where  no  such  detonation  solution  occurs. 
Additional  calculations  will  indicate  the  sensitivity  of  the  initial 
porosity,  <J>  ,  chemical  energy,  E  ,  burni ng  rate  index,  n,  and  ignition 

O  wiitM 


51 


energy. 


on  the  run-up  length  to  detonation. 


4. 2  Calculations 

Figures  17a  and  17b  present  a  case  where  a  transition  from  deflagra¬ 
tion  to  detonation  has  occurred.  For  this  example,  the  burning  rate  in¬ 
dex  was  n  -  1.0,  the  chemical  energy  E^^  =  5.48  MJ/kg  and  the  particle 

radius  r  =  127  um  (.005  in). 

P 

The  pressure-distance  profiles  at  five  separate  times  after  ignition 
of  the  propellant  at  x  =  0  are  shown  in  Figure  17a.  Examining  the  pro¬ 
file  for  t  =  50  ysec,  one  can  observe  that  the  profile  is  characterized 
by  a  shock  front  at  x  =  23  cm  followed  by  a  smooth  expansion  back  to  the 
wall  at  x  =  0  cm.  The  pressure  in  front  of  the  shock  is  at  atmospheric 
conditions  and,  therefore,  negligible  with  respect  to  the  pressure  behind 
the  shock. 

The  ignition  locus  plot  for  this  particular  case  is  shown  in  Figure 
17b.  Here,  the  ignition  front  moves  through  the  bed  at.  a  low  subsonic 


velocity  for  the  initial  ten  microseconds  and  then  accelerates  to  reach 
a  steady  state  velocity  of  =7.2  mm/ysec.  This  occurs  within  12  to  15 

cm  from  the  ionitprl  f»nd  .  Note  that  th^sp  df*tnnat  i  on  solutions  for  P 

and  are  in  fair  agreement  with  the  approximations  made  in  Ch ?  2 
(Eqs.  30  and  31)  for  the  given  input  »  V  Yj)- 

Obviously,  hydrodynamic  steady  state  analysis  (like  Eqs.  30  and  31) 
cannot  guarantee  that  a  transition  from  deflagration  to  detonation  will 


occur.  However,  it  seems  that  "critical"  values  of  porosity  (related  to 
gas  confinement),  gas  generation  rates,  and  chemical  energy  will  provide 
for  a  DDT.  For  example,  when  the  burning  rate  index  was  lowered  to 
n  =  0.8,  as  shown  in  Figures  18a  and  18b,  the  steep  pressure  front  associated 
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with  a  detonation  did  not  develop.  Correspond! ngly?  no  detonation  speed 
was  predicted.  In  this  example  the  peak  pressure  in  the  bed  never  ex¬ 
ceeded  5  GPa,  no  shock  was  predicted,  and  only  a  steady  convection-dr’ ven 
front  of  2.2  mm/usec  occurred  at  100  ysec  after  ignition  of  \  =  0.  It  is 
not  clear  at  this  time  whether  or  not  this  can  be  defined  as  a  "low-ve¬ 


locity''  detonation. 

According  to  Equation  30  in  Chapter  2,  for  the  case  where  a  transi¬ 
tion  to  detonation  actually  occurs,  the  steady  state  detonation  pressure, 

P  ,  should  increase  linearly  with  the  chemical  energy,  L  .  Also,  Equa - 

HIlM 

tion  31  states  that  the  detonation  velocity,  D^.  ,  is  a  function  of  the 
square-root  of  the  chemical  energy.  Figures  19a  and  19b  present  the  re¬ 
sults  of  a  calculation  where  all  parameters  were  identical  to  those  used 
to  give  the  DDT  results  of  Figures  17a  and  17b,  except  *  6.85  MJ/kg, 

an  increase  of  25%.  The  steady  state  shock  pressure  predicted  for  this 
case,  shown  in  Figure  19a,  was  21  GPa.  This  represents  a  (21/16.4  =  1.28) 
28%  increase  in  pressure  over  the  first  case.  The  predicted  detonation 
speed  (Fig.  19b),  calculated  from  the  slope  of  the  x-t  diagram,  was  8.70 


nun/usec.  According  to  Equation  31,  the  ratio  of  for  the  case  given  in 
Figure  13b  to  that  presented  in  Figure  17b  should  be 


A6.8S75.48)  -  1.32 

This  is  approximately  the  increase  predicted. 

Table  2  summarizes  the  detonation  pressure,  P  ,  and  the  detonation 
velocity,  D.T,  (which  were  the  end  results  of  the  TUT  calculations)  all 
as  a  function  of  the  propellant  chemica1  energy.  These  predicted  condi¬ 
tions  are  compared  with  the  approximate  steady  state  hydrodynamic  solu¬ 
tions  discussed  earlier,  i.e.,  Eqs.  30  and  31.  The  excellent  agreement 


'Location  (cm)  %  Press 


.g,  18a.  Pressure  History  During  Deflagration  with  no  Transition 
(n  =  0.8)  . 


Fig,  18b.  Ignition  Front  Locus  with  no  Detonation  [n  =  0.3). 


(cm } 


Fig.  19a.  Pressure  History  Leading  to  Detonation  with  25%  Increase 
THEM 


’i:  u-onijiiie  lo  case  shown  Vii  i  ig.  17a) 
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of  tho  detonation  pressure  with  the  analytic  solution  should  be  noted. 

However,  the  predicted  value  for  the  detonation  velocity,  D^,  from 
the  hydrodynamic  solution  (Eq.  31)  shows  to  be  slightly  less  than  the 
value  predicted  by  the  code  for  all  cases.  Although  one  cannot  judge 
which  of  the  two  values,  detonation  velocity  or  detonation  pressure,  is 
more  accurate,  the  percent  increase  in  the  detonation  velocity  as  the 
chemical  energy  is  increased  compares  favorably  with  the  hydrodynamic 
solution. 

4 . 3  DDT  Run-Up  Length 

The  run-up  length  to  detonation  is  defined  in  this  report  to  be  the 
distance  from  the  closed  end  where  the  bed  is  ignited  to  the  location 
where  both  the  peak  pressure  and  the  detonation  speed  are  constant,  i.e. 
the  equilibrium  steady  state  solution. 

Figures  20a  and  20b  plot  the  predicted  run-up  1 
as  a  function  of  the  burning  rate  pressure- index ,  n  (Fig.  20a),  and  the 
initial  bed  porosity,  4>Q(Figure  20b).  The  "no-solution"  boundary  indi¬ 
cates  that,  with  the  integration  scheme  used,  the  mesh  size  would  have 
to  be  drastically  reduced  (thereby  decreasing  the  integration  time  incre 
ment)  in  order  to  obtain  a  stable  solution.  This  is  a  costly  exercise, 
but  future  work  is  planned  to  increase  the  solution  regions.  Figure  20a 
clearly  indicates  that,  as  expected,  one  cannot  achieve  detonation  if 
che  burning  rate  during  the  deflagration  phase  is  not  sufficiently  large 
(see  "no  transition"  boundary).  For  the  solids  loading  considered 
(1-  <J)0  -  0.70),  it  would  appear  that  a  minimum  DDT  run  up  length  is  5  cm 
for  particles  of  250  pm  in  diameter. 

Figure  20b  begins  to  resemble  the  required  "U-shaped"  curve  of  l  (.  , 
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TABLE  2 

Comparison  of  the  Predicted  Detonation  State  with 
an  Approximate  Hydrodynamic  Solution 


echem 

V’CJ  (predict  1) 

PCJ  (Eq.  30) * 

D  (predicted) 

L  J 

DCJ(Eq.  31)** 

£CJ 

4.11  MJ/kg 

14.0  GPa 

12.0  GPa 

- 

6.48  -22- 
psec 

r  ,  _  mm 

5 .  Iz  - 

psec 

89  mm 

5.4S 

16.6 

16.0/ 

7.25 

5.92 

51 

6.85 

21.8 

20.08 

8.24 

6.62 

44 

8.22 

24.6 

warn 

^  -  1 7  |  /  .  2  j 

i 

7  O 

•J  o 

10.96 

30.0 

33.2 

10.58  !  3.37 

! 

■ 

n 

u 
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versus  initial  porosity.  DDT  experiments  by  Korotkov  et.  al.  [21]  show 
a  similar  behavior  (see  Fig.  21).  One  would  expect  that  for  a  relatively 
porous  bed,  4>q  >  0.60,  no  transition  will  occur  since  local  pressure  con¬ 
finement  is  limited.  If  the  porosity  is  too  small  no  gas  penetration 
for  the  accelerating  deflagration  wave  will  occur.  The  net  result  is  a 
porosity  where  a  minimum  run-up  length  occurs. 

Figure  22  presents  a  study  on  the  effect  of  the  chemical  energy,  Erucu, 

LilLM 

on  the  run-up  length,  ,  the  values  which  were  presented  in  'Iable  2. 

As  expected,  the  run-up  length  to  detonation  increases  as  the  amount  of 
chemical  enerQy  decreases.  Again,  as  expected  there  is  a  minimum  value 
where  no  transition  occurs,  =  3.0  MJ/kg.  For  this  case  it  does  ap¬ 

pear  that  a  constant  but  small  run-up  distance  is  still  required  as  the 
chemical  energy  increases  beyond  the  values  studied  here. 

4 . 4  Detonation  Reaction  Zone 

Gne  measure  of  the  fact  that  detonation  occurs  is  a  plot  of  the 
reaction  zone  width  versus  the  locus  of  the  ignition  front.  This  is  pre¬ 
sented  in  Figure  23,  where  the  reaction  zone  is  defined  as  the  region 
where  particles  are  ignited  and  generating  gas  (i.e.,  are  not  burned  out). 
As  shown,  the  zone  initially  increases  during  the  deflagration  phase  as 
the  convective  heat  transfer  provides  energy  to  ignite  more  and  more  of 
the  bed.  The  zone  then  collapses  to  a  thin  (constant)  width  as  the  sur¬ 
rounding  high  gas  presssure  causes  the  part'cles  to  burn  out  rapidly.  A 
steady  react! on  zone  thickness  of  approximately  y  mm  is  predicted.  How¬ 
ever,  most  of  the  gas  is  generated  in  a  small  region  immediately  behind 
the  ignition  front  where  the  particles  are  still  relatively  large,  .iote 
that  in  all  the  cases  reported  here,  the  initial  particle  diameter  was 
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2.S0  pm.  Obviously,  smaller  particles  will  provide  for  a  thinner  detona¬ 
tion  reaction  zone.  Figure  24  presents  a  porosity-distance  profile  for 
a  case  where  DDT  occurred.  Here,  a  porosity  of  <J>  =  Q.95  is  for  all  in¬ 
tents  and  purposes  the  condition  when  the  propellant  is  burned  out,  i.e., 
no  generation  of  hot  gases. 

4.5  Comments  and  Interpretations 

The  results  shown  in  Figures  17-24  have  clearly  indicated  that,  as 
expected,  high  solids  loading  of  relatively  small  particle  size  ener¬ 
getic  propellant  with  perfect  confinement  will  transit  into  a  detonation 
in  space  domains  of  several  centimeters.  One  would  expect  that  propel¬ 
lant  properties  and  packing  configurations  have  limits  where  detonation 
cannot  occur  and  this  was  clearly  shown  in  some  of  the  figures  which  show 
the  run-up  length  versus  property  parameters.  In  conclusion,  it  will  be 
useful  to  review  these  studies  to  determine  the  properties  and  configura¬ 
tions  which  minimize  a  DDT  hazard. 

For  example,  for  a  fixed  chemical  energy,  E^^  =  5.43  MJ/ kg,  a  fixed 

ignition  energy,  Ah.  =  9.0  KJ/kg,  a  fixed  particle  radius,  r  ■  127  pm 

^  k  ^ 

tO. 0,05  ini  and  a  fixed  solids  loadinv.  -fl-rtil  -  0.7.  the  burn! nv  rate  index. 

,  must  be  larger  titan  n  =*  G.St>  if  DDT  is  to  occur.  This  was  shown  in  Figure 
20a.  Of  course,  had  the  burning  rate  coefficient,  b,  been  a  different  num¬ 
ber,  this  exponent  may  have  been  different.  A  general  statement  >.n  then 
be  made  that  the  burning  rate,  and  hence  rate  of  gas  generation,  should  be 
kept  as  low  as  possible  to  minimize  DDT.  Conversely,  the  higher  the  burn¬ 
ing  rate,  the  better  the  chance  of  a  transition  to  detonation  occurring. 

Figure  20b  showed  the  run-up  length  to  detonation  versus  initial 
porosity,  <(,  .  As  was  stated  in  the  previous  section,  for  all  parameters 


.  24.  Porosity-Distance  Profile  for  a  Case  where  DDT  Occurred. 
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equal,  there  is  a  maximum  initial  porosity  where  no  transition  occurs. 

Since  a  randomly  packed  bed  of  unisized  particles  has  an  initial  porosity 
of  the  order  <J>o  =  0.4,  this  figure  indicates  that  this  is  in  the  region 
where  DDT  potential  is  at  a  maximum.  As  one  reduces  the  initial  porosity 
(that  is  increases  the  solids  loading)  the  run-up  length  to  detonation 
slightly  decreases  until,  as  experimental  work  of  Bernecker  and  Price  [9] 
has  indicated,  there  is  a  mimimum  initiax  porosity  where  DDT  cannot  occur. 
Since  at  these  initial  solids  loadings  multi-sized  particles  and  mechanical 
packing  are  required,  these  lo  lings  are  not  of  interest,  whereas  the  ran¬ 
domly  packed  loadings  are  of  interest. 

A  similar  comparison  of  the  run-up  length  to  detonation  was  shown  in 
Fig.  22  where  the  chemical  energy  content  was  varied.  Recall  that  a  chemi¬ 
cal  enerflv.  =  4.  IB  MJ/ke  11000  cal/el .  can  be  considered  an  ener- 

getic  propellant  material.  Based  on  the  results  shown  in  Fig.  22,  one 
can  state  that  less  energetic  material  than  this  has  little  chance  of  en¬ 
countering  a  DDT.  Doubling  the  chemical  energy  from  E^^  =  4.18  MJ/kg  to 
ECHF.M  =  8-36  MJAg,  represents  a  reduction  in  the  run-up  length  of  only 
about  one  half.  To  summarize,  high  energy  propellant  of  the  nitramir.e 
family  where  4.18  MJ/kg,  definitely  fail  within  the  regime  of  a 

DDT  hazard  if  properly  confined. 

The  final  comparison  of  this  type  is  shown  in  Figure  25.  This  shows 
the  z'un-up  length  to  detonation  versus  the  particle  radius.  As  one  would 
expect,  there  is  a  maximum  particle  radius  (i.e.,  surface-to-volume  ratio) 
where  transition  ls  not  predicted  to  occur.  This  is  indicated  by  the  "no- 
transition"  boundary  on  that  figure.  The  figure  also  indicates  that  as 
the  particles  get  smaller  in  size,  the  run-up  length  to  detonation  also 
decreases,  as  expected.  Reducing  the  initial  particle  size,  r  ,  to 
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values  less  than  25  pm  results  in  gas  generation  rates,  per  unit-volume, 

that  are  so  large  finer  grid  spacing  must  be  utilized  to  assure  stability. 

This  expensive  task  has  been  delayed  and  is  recommended  only  after  improved 

numerical  integration  schemes  have  been  developed. 

A  final  topic  studied  dealt  with  the  effect  of  ignition  temperature 

(or  more  appropriately  AE.  )  on  the  run-up  length  to  detonation.  For 

the  results  shown  in  this  chapter  a  nominal  value  of  T.  =  545°R  was  used. 

ign 

For  this  study  =  0.03,  rpQ  =  127  pm  and  n  =  1.0. 

Calculations  were  made  in  which  T.  varied  over  the  range  530°U  < 

ign 

T.  <  575°R.  For  a  constant  initial  bed  temperature  of  T  =  T  =  530°R, 
ign  -  p  g0  r0 

this  represents  a  range  of  ignition  energy  of  0.0  <  <_  25>, 

For  most  of  the  ignition  temperatures  tested  there  was  little  change  in 


the  steady  state  detonation  pressure  or  velocity.  Only  when 
approached  the  improbable  value  of  AE^  =  0-0  KJ/kg  did  the  values 
change  significantly. 
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APPENDIX  A 

JUMP  CONDITIONS  FOR  TWO-PHASE  REACTIVE  FLOW 


In  Section  1.3  (Jump  Conditions)  of  the  text  the  reader  was  provided 
with  a  review  of  the  jump  conditions  across  a  shock  discontinuity  for  one¬ 
dimensional,  one-phase  flow  with  heat  addition.  This  Appendix  will  out¬ 
line  the  development  of  the  jump  conditions  for  one-dimensional,  two-phase 
flow  with  heat  addition.  As  d  scussed  in  Section  2.3,  the  heat  addition 
for  two-phase  flow  comes  as  a  release  of  chemical  energy  from  the  ignited 
propellant  particles  to  the  surrounding  ga..c;.  For  a  detonation  this  is 
assumed  to  occur  in  an  infinitesimally  thin  reaction  zone. 

In  order  to  evaluate  the  jump  conditions  for  twe-phase  flow  across 
a  steady  state  combustion  wave,  the  conservation  equations  (Eqs.  22-27) 
are  first  written  as  conservation  equations  for  the  mixture. 

Mixture  Continuity: 

— -  [<f>p  u  ]  +  4 —  [(l-$)p  u  ]  =  0  (A.  1) 

ax  L  g  gJ  dx  p  pJ  J 

Mixture  Momentum: 


[$0  u2  +  4P  j  +  T-  [(!-<)>) P  +  (l-4>)  P  ] 
dx  g  gJ  dx  P  P  pJ 


=  0 


(A.  2) 


Mixture  Ej  ergy: 


4-  [<J>P  U  E  +  <J>U  P  ]  +  ~r~  [  ( 1  -d>)  p  U  E 
dx  L^g  g  gT  "  g  S  dx  11  p  Pl 


+  (l-4>Iu  P  ]  =  0 
P  P 


(A.  3) 


As  in  Secti,.->  1.3.  the  unreacted  or  "cold"  end  is  denoted  by  subscript  A 
and  the  products  or  "hot"  end  is  denoted  by  subscript  B  in  the  analysis 
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that  follows. 

The  mixture  conservation  equations  (Eqs.  A.1-A.3)  are  now  integrated 
froi  the  "cold"  end  to  the  "hot"  end  of  the  bed  with  the  boundary  conditions 
<KB)  =  1.0 

VA) "  V*1  ■  VA 

V“>  •  VB 

PS(B)  =  pb 

Integrated  Mixture  Continuity: 


PBVB  -  VgA  VA  +  C1^A)ppAVA 


Integrated  Mixture  Momentum: 


p»  -  *  p„  J 

-A 


(A.  4) 


*  "  V">pava  ♦  V 


Intcg’. r  •  t Mixture  Lb.eigy: 

gT 


,pv  vl  *  V«  “  Kvv! 


(A.  5) 


(l-itip  v(L:v-  ;  i  * 
h  'j  jA 


i’jVj  "  [Cl-tn’pV 


(A.  6) 


At  first  glance.  Equation  A. 4  -'o  h..f  ;eem  to  :  -.asily  evaluated  know¬ 
ing  the  conditions  at  both  ends  of  the  flow.  However,  as  was  pointed  out 
by  Kuo  and  Suinmerf  ield  [23],  they  are  not  simple  ilgcbTaic  equations  un¬ 
til  P  ,  the  stress  transmitted  through  the  parti-;'-  phase  at  location  A. 

Fa 

is  evaluated.  This  was  stiown  to  be  obtainable  by  •  .tegratlng  the  momentum 
equation  from  tin.  cold  end.  A,  to  the  location  wi  >» .  e  the  particles  are  no 
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longer  in  contact,  denoted  by  subscript.  C.  From  this  you  obtain 

fC  C 

-P  (1-$.)  -  -  V  dx  +  ”  u  dx 

PA  A  P 


tC  rC 

V  dx  +  ”  u  dx 

P 

A  A 


d  u*  (1-40 

V2(l-<j>) 

Lp  p  Jc  L 

p  -1a 

It  should  be  noted  that  if  the  flow  situation  being  modeled  consists 

of  an  infinttly  long  bed,  then  the  particle  stress  ac  location  A,  P  , 

PA 

will  also  be  equal  to  zero.  This  is  true  since  pressure  disturbances  ini¬ 
tiated  at  the  combustion  zone  (i.e.,  compression  of  the  propellant  par¬ 
ticles)  are  propagated  through  the  medium  at  the  speed  of  sound  of  the 
solid,  and  would  therefore  take  an  infinite  amount  of  time  to  reach  loca¬ 
tion  A. 

Since  all  experimental  work  consists  of  finite  length  beds  (the  nu¬ 
merical  model  studied  in  this  report  was  L  =  25  cm) >  a  comment  on  the 
integration  limits  of  Equation  A. 7  is  appropriate.  If  ,  te  region  where 
particle  drag  is  of  significance  extends  from  the  combustion  wave  (C)  to 
the  wall  condition  (A),  then  the  integrals  in  Equation  A. 7  must  be  evalu¬ 
ated  after  each  time  increment  and  hence  the  particle  stress  at  the  wall, 

P„  ,  is  changing.  However,  if  the  region  of  significant  particle  drag 
^A 

is  assumed  to  be  only  that  z  >ne  immediately  ahead  of  the  combustion  zone, 

then  the  integrals  in  Eq.  A. 7  are  always  constant  and,  therefore,  P  is 

PA 

also  constant. 

A  final  comment  on  the  particle  stress  is  that  if  the  time  to  detona¬ 
tion  is  less  than  the  time  it  takes  a  pressure  disturbance  to  travel 
through  the  solid  matrix  to  the  rear  wall,  then  the  rear  wall  has  not 
felt  the  pressure  disturbance  and  hence,  P  -  0. 


0. 
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APPENDIX  B 


CONSTITUTIVE  RELATIONS 


As  discussed  in  Chapter  2,  a  system  of  nine  independent  equations  is 

necessary  in  order  to  solve  for  the  nine  unknown  variables  describing  two- 

phase  flow:  p,p,u,u,E,E,P,P  and  4>. 

F  g  P  g  P  g  P  g  P 

Cf  the  three  additional  equations  necessary  for  the  solution,  one  is 
the  nonideal  equation  of  state  for  the  gas  phase  as  described  in  the  text. 


P 

_JL_ 

RT 

g 


1  +  bp 

g 


c  (bpg) 2  +  .  .  . 


(B.l) 


Values  for  the  constants  b  and  c  are  discussed  in  Appendix  C. 

The  second  represents  an  equation  of  state  which  relates  the  solid 
density  to  the  stress  on  the  particle.  (See  Ref.  3).  This  is  also  men¬ 
tioned  in  the  text. 


p  3  K 

3  -  ((-2-)  -  i)  ~ 

p  “p  J  J  3 


(B.2) 


Finally,  the  third  additional  constraint  is  a  relation  for  the  par¬ 
ticle  phase  stress,  P  ,  as  a  function  of  the  solids  l.-ading  and  the  material 

bulk  modulus,  K  (also  see  Ref.  3J . 

o  J 

As  discussed  in  References  2  and  3,  one  must  also  specif>  functional 
relations  for  the  following: 


a)  1’  =  mass  generation  rate  per  unit  ’'olume 
T  =  — •  (l-it)  p  r 


(B.3J 


Here,  r 

p 


the  instantaneous  panicle  radius  and 


is  t he  surface  burning 


75 


rate  specified  as  a  function  of  pressure  (and  possible  partible  temperature ) , 
For  all  cases  run  in  this  study. 


r  =  b" 


CE.4) 


where  r  has  units  of  (in/sec),  P  (psi)  and  _b  is  of  the  order  (1  x  30'') 
in/sec  . 

(psi) 

b)  V  =  interphase  (gas-particle)  viscous  force  (as  discussed  ir.  the 
text) . 


where 


V  =  (u  -  u  )  f 

4r  2  8  P  P8 

P 


f  =  5.C6  x  105  r  Re/u  Z  (see  Eq.  37) 
P2  P  g 


(3.5) 


c)  Q  =  interphase  heat  transfer  rate 


Q  =  —  (1-40  h  (T  -  T  ) 
rp  PR  g  P 

In  the  analysis  carried  out  here,  the  heat  transfer  coefficient  was 


(B.6) 


h  =  0.65  [Re)U,/  (Pr) U 

P  P 


(B .  7) 


where  k^  is  the  thermal  conductivity,  Re  is  the  Reynolds  number  and  Pr  is 


the  Prandtl  number . 
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APPENDIX  j 


NON I DEAL  EQUATIONS  OF  STATE 


As  described  bv  Jac  obs  |  Uj ,  the  nonideal  equation  of  state  is  ini¬ 
tially  written  ns  a  polyru.ini a i  expansion  in  gas  density. 


Pv 

RT 


1  +•  x 


a 

ex"  -*•  d> 


CC.l) 


Here,  x  =  bp  -  — . 

v 

If  Cv>  the  specific  heat  at  a  constant  volume,  is  assumed  constant 
over  the  varying  temperatures,  as  was  in  the  model,  Eq.  C.l  can  be 
written  as: 


Pv  =  (^-1)  E  f(v)  (C.2, 

C 

In  Eq.  C.2,  V.  is  the  specific  heat  ratio  for  the  ideal  gas  limit,  y.  --  , 

1  v 

and  f(v)  is  the  gas  density  polynomial  expansion. 


, . ,  ■  ,  b  ,b .  “ 

i  (v  i  -  1  +  —  +  c  (—)  *• 

'  v  v 


(C 


For  simplicity  t-ie  higher  order  terms  are  dropped. 

The  local  speed  of  sound,  a,  in  the  detonation  state  is  defined  as: 


2 

a 


Here,  y  is  the  effective  specific  J. sat  ra-'o  in  t lie  detonation  state. 
For  detonating  •••  j  lusise:-  ■  fvp1.  ■  ’  1\  we<m  >  v..  to  •  h  ee  .  The 

equality  in  Equation  .  4  ■  ■■  no-.-  s<  t  :  r<  ,<■ 
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Y  =  £  (—)  =  £  (IB.)  (l£)  -  i i  (IB.)  r.  .L_) 

Y  P  C3pJ  P  '■3v'  ’■Sp'1  P  l3vJ  *■  ' > 1 


s  p 


o: 


“V  ,3P, 

v  =  P  (9vj 


(C.5) 


Since  P  --  P  (fc ,  v)  , 


d!  *  a!.:  -  (~)  dv 


(C.6) 


and 


dP  J»j_ ,  di 

^d  v  3E'  d\ 


i.\  ,dv, 

'dv  _  ‘-dv 


(C.7) 


Substituting  Eq  .  C.‘:  an--  tne  constant  entropy  thermodynamic  relation 


P  =  ■ 


p-.o  F.q 


-  j 

:T 


':*P) 


(c.») 


f:  ''  -t,  t:  -  ; 


•  T  Li  ..r-  * 

I  of  eq  S  ?  a. 
•7  )  ■ 


M/af  •••**•  a  rm<»ar  i  n  a  i  ri  Fn  C  -  ^ 


v  ..s  taken  holding  t  constant, 
f’(v) 


•3v 


(C.  10) 


and  ‘,  conu  « rt  a  *.  i. 


r '  o  t:  is  taken.  adding  v 


l-L 


f  * 

iM  — .  —  • 


(c; .  ]  J.  : 


m.  •  si! 


i*.  n  j  j* 


,  i>  r.».  •»*  .a> 
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the  form: 


y 


(Yi-13  f(v) 


1Y.-1)  E  f(v) 


CY.-l)  E  f '  (v)  v 
Pv 


But,  since  from  Eq.  C.2 


CC. 12) 


Pv  1 
(Yj-DE  "  f(v) 

Eq.  C.12  is  now: 

Y  =  (Y-l)  f(v)  -  1  -  v  f-^p-  (C.  13 j 

Knowing  approximately  what  value  of  y  an  explosive  exhibits,  values  of 
the  constant  b  and  c  can  be  fit  to  equate  Eq.  C.13. 
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APPENDIX  D 


FINITE  DIFFERENCING  TECHNIQUES 


This  section  will  give  the  reader  a  more  detailed  description  of  the 
finite  differencing  scheme  referred  to  in  Chapter  3,  the  Lax-Wendroff 
center  differencing  technique. 

In  order  to  solve  the  system  of  six  hyperbolic  partial  differential 
equations  describing  two-phase  reactive  flow,  the  Lax-Wendroff  scheme  pre¬ 
sented  in  Reference  19  for  one-phase  flow  was  modified  to  model  two -phase 
flow.  This  modification  was  necessary  because  of  the  interactive  source- 
sink  terms  involved  in  the  governing  equations. 

The  six  conservation  equations  described  in  Chapter  2  (Eqs.  22  to  21) 
can  be  written  in  vector  form  as: 


30  t  3F 
9t  *  dx 


(D.l) 


In  this  notation  U,  F  and  C  are  the  vectors 


P1 

P2 

°2Up 

P1UH 

(pi%  *  *  y 

P,u 

2  P 

F  = 

tP2Up  *  Cl~<}>)Fp) 

PlEgT 

(p.u  E  +  u  PI 

1  nT  g  g 

P2EgT 

P2UPLPt  ♦  (1-<J>)  Up  P  ) 
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C  = 


r 

-r 


r  u  -  v 

p 

-r  u  +  v 

p 


r‘E»EM  *  -$>  '  ®  S  '  « 


r^ECHEM  ^  ^  13  up  *■  Q 


The  Lax-Wendrof f  solution,  as  described  by  Richtmeyer  and  Morton  [22] 
begins  as  a  Taylor's  series  expansion  accurate  to  the  second-order. 

\  2.  t 


,.t+At  „t  A4.  [Sul*  1  .  2  [3  u 

U  =  U  +  At  -s—  +  -=■  At  —rr 

x  X  at  I  2  U  2 

*-  -ix  -ot  - 


(I).  2) 


In  order  to  make  use  of  Equation  V.2,  the  time  derivatives  must  be  replaced 
by  derivatives  in  space. 

3Fi 

By  using  Equation  D.l  and  introducing  the  matrix  A,  where  A  .  =  ~ 


3  II  3_ 
2  =  3t 

at 


1 - 1 

t 

X  |*TJ 

3C 

’3Fl 

3  C 

3__ 

r 

3U  ' 

3C 

3_  [  3F 

=  3t 

3x 

_3tj 

3  3c  " 

3x 

la 

at  _ 

at  ‘ 

3x  3x 

.  m 

CD. 3) 


Replacing  x-derivat ives  with  finite  difference  quotients  (i.e., 

3F/3x  =  (Fx+Ax  -  Ax3/2Ax),  the  Lax-Wendroff  solution  can  now  be  written 


,t»-At 


u*  -  -A-i 

x  2Ax 


Fl  .  -  F*  .  j  +  AfC 

x+Ax  x-Axj 


[At]2 

'  t 

rt  ,.t  1 

.  t 

r  t  t  V 

[AxJ 

x+Ax 

L  2 

F  .  -  F 

x+Ax  xj 

-  A 

|  x- Ax 

2 

F  -  F  . 
x  x-Ax 

+ j  At *C  (D.4) 


r 
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A*"  .  denotes  A  [x  .  +4  ll*]  ■  For  A  a  constant  matrix  the  Lax-Wendroff 

x+Ax  2  x+Ax  2  xJ 

£. 

solution  takes  the  form 


CD. 5) 


Richtineyer  [22]  developed  a  variation  of  the  Lax-Wendroff  solution  al¬ 
so  of  second-order  accuracy.  This  method  involves  two  steps  and  eliminates 
the  use  of  the  A  matrix  found  in  the  Taylor's  series  expansion  solution. 

Richtmeyer's  first  step  involves  computing  intermediate  values  in 
time  and  x- location  of  U. 


„  At 
t+  -s- 

J  * 

X+  Ax 

2 


■  1  K'x.ix  *  <  ]  *  m  [  rx.Ax  -  Fx  ]  '  r  •  cx 


The  first  step  (.Eq.  U.t>j  is  of  first-order  accuracy.  However,  the  final 
step  (Eq.  D.7)  is  accurate  to  the  second-order  since  r  quantities  of 
8  (Ax)  are  differenced  over  Ax.  For  the  linear  case  where  F(U)  =  AU 
Equations  D.6  and  D.7  combine  to  yield  the  Lax-Wendroff  solution,  Equa¬ 
tion  D.5. 

A  two-dimensional  (x-t)  grid  is  helpful  in  understanding  how  the 
Lax-Wendroff  differencing  schome  can  be  applied  to  the  numerical  solution 
■>\  the  conservation  equations  on  a  digital  computer.  In  Figure  D. 1  a  known 


S2 


value  of  the  vector  U  (mass,  momentum  and  energy  per  unit  volume)  at  a 

time  t  =  t  and  located  at  x  *  xq  is  indicated  on  the  grid  by  a  solid 

dot  at  (xq,  tQ) .  Locations  where  the  vector  U  is  yet  to  be  solved  at 

are  indicated  by  open  circles.  Figures  D.2  and  D.3  illustrate  the 

Richtmeyer  two-step  routine  using  the  x-t  grid.  For  instance,  in  Figure 

D.2  the  vector  U  at  (x  -  .  t  *  4^-1  is  solved  for  by  using  known  values 

O  2  O  2  ' 

at  two  locations  (x  ,  t  )  and  (x  -  Ax,  t  ) . 

v  o’  o'  v  o  o' 

In  a  later  paper  by  Rubin  and  Berstein  [19]  it  was  shown  for  cne- 

phase  flow  that  a  modified  version  of  Richtmeyer 1 s  2-step  method  was 

more  stable  in  solving  the  system  of  equations.  The  use  of  half-time 

steps  was  eliminated  by  averaging  F  differences  at  present  time,  tQ,  and 

future  time,  t  +  At.  This  is  written  as 
o 

Step  1 


,.t+At  _  x+Ax  *  x  At  ,t  „t] 

Ux+Ax  "2  +  Ax  x»Ax  *  rxJ 

2 


(.b.S) 


Step  2 


U 


t+At 

x 


At 

+  2Ax 


■  1 

F 

t  1  FC*At 

(,t*At.  ‘ 

2 

x*  Ax 

x- Ax  1  x+ Ax 

‘  x-Ax 

-1  T 

2  - 

-•9) 


This  method  is  illustrated  in  Figures  D  4  and  D. 5  on  the  following  page. 

The  following  is  a  summary  of  the  schemes  described  in  this  Appendix. 
Each  will  be  applied  to  the  gas  continuity  equation  of  l-D,  two-phase  flow. 
Here,  T  represents  the  gas  being  generated  from  the  solid  phase. 

Gas  continuity  equation 


at 


apiuf, 

ax 


(u.  10) 
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